rm(list= ls())

library(tidyverse)
library(fixest)
library(fastDummies)
library(lfe)
library(broom)

## Define rank_fun

rank_fun <- function (x) {
  su = sort(unique(x))
  for (i in 1:length(su)) x[x == su[i]] = i
  return(x)
}

# NRW:
# left_total = SPD + Greens + Left
# totright2_nofdp = CDU/CSU + AfD 

nrw <- read_rds("data/nrw_elections_clean.rds")

nrw <- readRDS('data/state_elections_clean.RDS') %>% 
  filter(land == '05') %>%
  dplyr::select(1:15) %>% 
  mutate(year = lubridate::year(date)) %>%
  mutate(post = ifelse(year == 2017, 1, 0)) %>%
  mutate(totright2_nofdp = cdu_csu + ifelse(year == 2017, afd, 0)) %>% 
  mutate(period = rank_fun(year))

## Need to get treatment

control_counties <- readRDS('data/treated_by_county.rds') %>%
  filter(state == 'Nordrhein-Westfalen' & treated == 0) %>%
  pull(county.id)

nrw <- nrw %>%
  mutate(treated = ifelse(county %in% control_counties, 0, 1),
         treated_post = treated*post)

## Merge-in AA District IDs 

aadist <- readRDS('data/treated_by_county.rds') %>%
  dplyr::select(AA.district.id, county.id)

nrw <- nrw %>% left_join(aadist, by = c('county' = 'county.id'))

## Nr of pre and post treatment periods

id_treat <- nrw %>% 
  filter(post == 1) %>% 
  filter(period == min(period)) %>% 
  pull(period) %>% 
  unique()

## Create Period*Treated Indicator Variables 

nrw <- fastDummies::dummy_cols(nrw, select_columns = 'period') %>% 
  mutate_at(vars(contains('period_')), ~ .*treated) 

## Generate Model Formulas 

outcomevars <- c('left_total', 'totright2_nofdp')
periodvars <- names(nrw)[str_detect(names(nrw), 'period_')]

## Exclude last pre-treatment period 

periodvars <- periodvars[-which(periodvars == 'period_2')]

results <- lapply(outcomevars, function(dv){
  
  model_formula <- as.formula(paste0(dv, ' ~', paste0(periodvars, collapse = '+'), '| period + ags |0 | AA.district.id'))
  
  model_est <- felm(model_formula, data = nrw) %>%
    broom::tidy(conf.int = T) %>% 
    filter(str_detect(term, 'period')) %>% 
    mutate(outcome = paste(dv))
  
}) %>% reduce(rbind) 

## Rename the period variables

results <- results %>% 
  mutate(period = parse_number(term) - (as.numeric(id_treat) - 1))

##

results <- results %>% 
  mutate(outcome = recode(outcome, `left_total` = 'Total left vote',
                          `totright2_nofdp` = 'Total right vote'))


## Add zero 

zero_period <- data.frame(term = rep('period_2', 2),
                          estimate = rep(0, 2),
                          conf.low = rep(0, 2),
                          conf.high =rep(0, 2),
                          outcome = unique(results$outcome),
                          period = rep(0, 2))

results <- bind_rows(results, zero_period) %>%
  mutate(period = as.factor(period))


## Plot this

myLoc <- 
  (which(levels(factor(results$period)) == "1") +
     which(levels(factor(results$period)) == "0")) / 
  2


p1 <- ggplot(results, aes(x = factor(period), y = estimate, 
                          ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = 'dotted', color = 'grey40') +
  geom_vline(xintercept = myLoc, linetype = 'dotted', color = 'grey40') +
  geom_errorbar(width = 0) +
  geom_point(fill = 'white', shape = 21) +
  theme_bw() + 
  ylab('Coefficient Estimate\n(percentage points)') + 
  xlab('') + 
  facet_wrap(~ outcome) + 
  scale_color_grey(name = '') +
  scale_x_discrete(breaks = as.character(c(-1, 0, 1)),
                   labels = c('2010', '2012', '2017')) +
  ylim(c(-8, 6))
p1

